library(dplyr)
library(MASS)
library(candisc)
library(klaR)
library(FSA)
library(ROCR)

library("PerformanceAnalytics")
library("FactoMineR") 
library("factoextra")

Подготовка данных

Возьмем набор данных Wine quality dataset.

Этот набор данных включает данные о физико-химическом составе красных вин.

df_raw <- read.csv("winequality-white.csv", sep = ";")
head(df_raw)

Признаки:

Количество индивидов:

length(df_raw[, 1])
## [1] 4898

Посмотрим на количество классов.

count(df_raw, quality)

И на их соотношение:

perc_class <- count(df_raw, quality)
perc_class$n <- count(df_raw, quality)$n / sum(count(df_raw, quality)$n)
perc_class

Обозначим классы qualuity:

df <- df_raw %>% filter(quality > 4, quality < 8)
df_low <- df_raw %>% filter(quality < 5)
df_low$quality <- 5
df_high <- df_raw %>% filter(quality > 7)
df_high$quality <- 7
df <- rbind(df, df_low, df_high)
df$quality <- cut(df$quality, 3, labels=c('1', '2', '3'))

Проверка признаков

Посмотрим, есть ли признаки, которые стоит логарифмировать.

chart.Correlation(subset(df, select=-c(quality)), histogram=TRUE)

Прологарифмируем residual.sugar, volatile.acidity, citric.acid, chlorides, free.sulfur.dioxide, sulphates, density.

df_l <- df

df_l <- transform(df_l, residual.sugar=log(residual.sugar))
df_l <- transform(df_l, chlorides=log(chlorides))
df_l <- transform(df_l, citric.acid=log(citric.acid))
df_l <- transform(df_l, volatile.acidity=log(volatile.acidity))
df_l <- transform(df_l, free.sulfur.dioxide=log(free.sulfur.dioxide))
df_l <- transform(df_l, sulphates=log(sulphates))
df_l <- transform(df_l, density=log(density))

chart.Correlation(subset(df_l, select=-c(quality)), histogram=TRUE)

Есть 19 индивидов, у которых citric.acid был равен нулю. Удалим эти индивиды.

df_l <- df_l %>% filter(!is.infinite(citric.acid))

Train Test split

Разобьем выборку на \(train\) и \(test\).

smp_size <- floor(0.75 * nrow(df))

set.seed(123)
train_ind <- sample(seq_len(nrow(df_l)), size = smp_size)

train <- df_l[train_ind, ]
test <- df_l[-train_ind, ]

rownames(train) = 1:length(train[, 1])
rownames(test) = 1:length(test[, 3])

Посмотрим на баланс классов в \(train\) и \(test\).

count(train, quality)
count(test, quality)

И на их процентное соотношение:

perc_class <- count(train, quality)
perc_class$n <- count(train, quality)$n / sum(count(train, quality)$n)
perc_class
perc_class <- count(test, quality)
perc_class$n <- count(test, quality)$n / sum(count(test, quality)$n)
perc_class

MANOVA

Проверим гипотезу о равенстве мат. ожиданий среди групп quality.

train.manova <- manova(cbind(fixed.acidity, volatile.acidity, citric.acid, residual.sugar, chlorides, free.sulfur.dioxide, total.sulfur.dioxide, density, pH, sulphates, alcohol) ~ quality, data = train)

summary(train.manova, 'Wilks')
##             Df   Wilks approx F num Df den Df    Pr(>F)    
## quality      2 0.67848   71.217     22   7320 < 2.2e-16 ***
## Residuals 3670                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(train.manova, 'Roy')
##             Df     Roy approx F num Df den Df    Pr(>F)    
## quality      2 0.44341   147.58     11   3661 < 2.2e-16 ***
## Residuals 3670                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Однако перед тем, как интерпретировать результаты критериев, посмотрим еще раз на распределение прологарифмированных данных. Попробуем понять хотя бы приблизительно, равны ли ковариационные матрицы среди групп quality.

ggplot(df_l) + 
  geom_point(aes(x = pH, y = alcohol, group = quality, color = quality)) + 
  stat_ellipse(aes(x = pH, y = alcohol, group = quality, color = quality)) + 
  scale_color_manual(values = c("red", "purple", "black"))

ggplot(df_l) + 
  geom_point(aes(x = fixed.acidity, y = volatile.acidity, group = quality, color = quality)) + 
  stat_ellipse(aes(x = fixed.acidity, y = volatile.acidity, group = quality, color = quality)) + 
  scale_color_manual(values = c("red", "purple", "black"))

По формам эллипсов можно уже заметить, что ковариационные матрицы различны, откуда следует тот факт, что на критерии Wilk’s Lambda или Roy’s greatest root нельзя полагаться.

Более того, из определения модели \(LDA\) плотность в точке \(x\): \(p_i(x) = p(x|\xi = A_i) = \frac{1}{2\pi^{p/2}|\Sigma|^{1/2}}exp(-\frac{1}{2}(x-\mu_i)^\mathrm{T}\Sigma^{-1}(x-\mu_i))\) и классифицирующие функции \(f_i(x) = \pi_ip(x|\xi = A_i)\) предполагают равенство ковариационных матриц между классами. Хотя мы и проверяем “на глаз”, скорее всего нам больше подойдет \(QDA\), нежели \(LDA\).

LDA + error tables

Обучим модель на train и на нём же протестируем.

full.lda <- lda(x=train[, 1:11], grouping = train[, 12])
full.ldap <- predict(full.lda, train[, 1:11])
full.ldapc <- full.ldap$class
ct <- table(full.ldapc, train[, 12])

ct
##           
## full.ldapc    1    2    3
##          1  700  332   30
##          2  491 1141  446
##          3   21  200  312

Точность модели по классам:

diag(prop.table(ct, 2))
##         1         2         3 
## 0.5775578 0.6820084 0.3959391

Попробуем задать априорные вероятности \(\pi_i\) равными, а не на основе соотношения количества элементов в выборке.

full.lda <- lda(x=train[, 1:11], grouping = train[, 12], prior = c(1, 1, 1)/3)
full.ldap <- predict(full.lda, train[, 1:11])
full.ldapc <- full.ldap$class
ct <- table(full.ldapc, train[, 12])

ct
##           
## full.ldapc   1   2   3
##          1 857 502  85
##          2 257 599 154
##          3  98 572 549

Точность (accuracy) модели по классам:

diag(prop.table(ct, 2))
##         1         2         3 
## 0.7070957 0.3580395 0.6967005

Видим что точность для \(1\)го и \(3\)го класса сильно возрасла, в то время как для \(2\)го класса упала.

Посмотрим на результаты при обучении на train и проверке на test.

Обучим модель:

train.lda <- lda(train[, 1:11], train[, 12], prior = c(1, 1, 1)/3)
train.ldap <- predict(train.lda, test[, 1:11])
train.ldapc <- train.ldap$class
ct <- table(train.ldapc, test[,12])

ct
##            
## train.ldapc   1   2   3
##           1 291 160  27
##           2  89 168  56
##           3  37 189 189

Точность модели по классам:

diag(prop.table(ct, 2))
##         1         2         3 
## 0.6978417 0.3249516 0.6948529

Построим модель используя кросс-валидацию:

cv.lda <- lda(df_l[, 1:11], df_l[, 12], CV = TRUE, prior = c(1, 1, 1)/3)
cv.ldac <- cv.lda$class
ct <- table(cv.ldac, df_l[,12])
ct
##        
## cv.ldac    1    2    3
##       1 1170  670  116
##       2  336  757  217
##       3  123  763  727
diag(prop.table(ct, 2))
##         1         2         3 
## 0.7182320 0.3456621 0.6858491

Canonical Analysis

result.candisc <- candisc(train.manova)

ggplot(result.candisc$scores) + 
  geom_point(aes(x = Can1, y = Can2, color = quality)) +
  stat_ellipse(aes(x = Can1, y = Can2, color = quality)) +
  geom_segment(as.data.frame(result.candisc$structure), 
               mapping = aes(x = 0, y = 0, xend = 6.5 * Can1, yend = 6.5 * Can2), 
               arrow = arrow(angle = 10, type = "closed"), size = 0.5)

Изображение разделяющих плоскостей в разных проекциях

plot(partimat(quality ~ ., train, method="lda"), xlim = c(-10, 10),
     ylim = c(-10, 10));

Можно заметить, что в большинстве графиков один из классов (фиолетовый) отделяется от остальных, мягко говоря, плохо. К тому же, индивиды сильно перемешаны.

ROC-AUC

Посмотрим на графики \(ROG\) попарно для каждых комбинаций \(2\)х классов.

Сначала возьмем класс \(1\) и \(2\).

# Подготовка данных
two.df = Subset(df_l, quality != 3)
two.lda <- lda(two.df[, 1:11], grouping = two.df[, 12], prior=c(1,1)/2)
two.ldap <- predict(two.lda, two.df[, 1:11])

# Переводит данные в стандартизированный формат. Первый параметр - предсказания, второй - настоящая классификация.
preds <- prediction(two.ldap$posterior[,2], two.df[, 12])

# Строим ROC-кривую и рисуем ее
perf <- performance(preds, "tpr", "fpr")
plot(perf, colorize = TRUE)

#Считаем AUC
AUC.ROCR <- performance(preds, "auc")

# Добавляем прямую
abline(a = 0, b = 1)

# Добавляем значение AUC на график
text(x = .25, y = .65, paste("AUC = ", round(AUC.ROCR@y.values[[1]],5)))

Теперь возьмем класс \(2\) и \(3\).

# Подготовка данных
two.df = Subset(df_l, quality != 1)
two.lda <- lda(two.df[, 1:11], grouping = two.df[, 12], prior=c(1,1)/2)
two.ldap <- predict(two.lda, two.df[, 1:11])

# Переводит данные в стандартизированный формат. Первый параметр - предсказания, второй - настоящая классификация.
preds <- prediction(two.ldap$posterior[,2], two.df[, 12])

# Строим ROC-кривую и рисуем ее
perf <- performance(preds, "tpr", "fpr")
plot(perf, colorize = TRUE)

#Считаем AUC
AUC.ROCR <- performance(preds, "auc")

# Добавляем прямую
abline(a = 0, b = 1)

# Добавляем значение AUC на график
text(x = .25, y = .65, paste("AUC = ", round(AUC.ROCR@y.values[[1]],5)))

И последняя пара - классы \(1\) и \(3\).

# Подготовка данных
two.df = Subset(df_l, quality != 2)
two.lda <- lda(two.df[, 1:11], grouping = two.df[, 12], prior=c(1,1)/2)
two.ldap <- predict(two.lda, two.df[, 1:11])

# Переводит данные в стандартизированный формат. Первый параметр - предсказания, второй - настоящая классификация.
preds <- prediction(two.ldap$posterior[,2], two.df[, 12])

# Строим ROC-кривую и рисуем ее
perf <- performance(preds, "tpr", "fpr")
plot(perf, colorize = TRUE)

#Считаем AUC
AUC.ROCR <- performance(preds, "auc")

# Добавляем прямую
abline(a = 0, b = 1)

# Добавляем значение AUC на график
text(x = .25, y = .65, paste("AUC = ", round(AUC.ROCR@y.values[[1]],5)))

Как мы видили на графиках в разделах ранее, классы \(1\) и \(3\) являлись самыми разделимыми, что мы и наблюдаем на \(ROC-AUC\).

Но в целом, модель плохо справляется с поставленной задачей.

QDA

Попробуем взять модель \(QDA\).

Обучим ее:

full.qda <- qda(x=train[, 1:11], grouping = train[, 12])
full.qdap <- predict(full.qda, train[, 1:11])
full.qdapc <- full.qdap$class
ct <- table(full.qdapc, train[, 12])

ct
##           
## full.qdapc   1   2   3
##          1 756 376  44
##          2 391 863 251
##          3  65 434 493

Точность модели по классам:

diag(prop.table(ct, 2))
##         1         2         3 
## 0.6237624 0.5158398 0.6256345

Изначально точность \(3\)го класса получилась выше по сравнению с \(LDA\).

Попробуем задать равные априорные вероятности \(\pi_i\).

full.qda <- qda(x=df_l[, 1:11], grouping = df_l[, 12], prior = c(1, 1, 1)/3)
full.qdap <- predict(full.qda, df_l[, 1:11])
full.qdapc <- full.qdap$class
ct <- table(full.qdapc, df_l[, 12])

ct
##           
## full.qdapc    1    2    3
##          1 1130  635   70
##          2  282  579  162
##          3  217  976  828

Точность модели по классам:

diag(prop.table(ct, 2))
##         1         2         3 
## 0.6936771 0.2643836 0.7811321

Видим что точность для \(3\)го класса сильно возрасла, в то время как для \(2\)го класса упала.

Посмотрим на результаты при обучении на train и проверке на test.

Обучим модель:

train.qda <- qda(train[, 1:11], train[, 12], prior = c(1, 1, 1)/3)
train.qdap <- predict(train.qda, test[, 1:11])
train.qdapc <- train.qdap$class
ct <- table(train.qdapc, test[,12])

ct
##            
## train.qdapc   1   2   3
##           1 280 143  16
##           2  81 133  39
##           3  56 241 217

Точность модели по классам:

diag(prop.table(ct, 2))
##         1         2         3 
## 0.6714628 0.2572534 0.7977941

Из-за маленького объема выборки, результаты не особо вселяют уверенность.

Построим модель и протестируем ее используя кросс-валидацию:

cv.qda <- qda(df_l[, 1:11], df_l[, 12], CV = TRUE, prior = c(1, 1, 1)/3)
cv.qdac <- cv.qda$class
ct <- table(cv.qdac, df_l[,12])
ct
##        
## cv.qdac    1    2    3
##       1 1117  642   71
##       2  293  561  173
##       3  219  986  816
diag(prop.table(ct, 2))
##         1         2         3 
## 0.6856967 0.2562814 0.7698113

Точность модели для \(2\)го класса оставляет желать лучшего.

plot(partimat(quality ~ ., df_l, method="qda"), xlim = c(-5, 5), ylim = c(-5, 5))

На практике получили точность у \(QDA\) слегка хуже, чем у \(LDA\), однако из-за того, что облака точек находятся близко друг к другу, ни одна из этих моделей не может более-менее точно классифицировать качество вин.

ROC-AUC QDA

Посмотрим на графики \(ROG\) попарно для каждых комбинаций \(2\)х классов. Будем смотреть сразу для модели с кросс-валидацией.

Сначала возьмем класс \(1\) и \(2\).

# Подготовка данных
two.df = Subset(df_l, quality != 3)
two.qda <- qda(two.df[, 1:11], grouping = two.df[, 12], prior=c(1,1)/2)
two.qdap <- predict(two.qda, two.df[, 1:11])

# Переводит данные в стандартизированный формат. Первый параметр - предсказания, второй - настоящая классификация.
preds <- prediction(two.qdap$posterior[,2], two.df[, 12])

# Строим ROC-кривую и рисуем ее
perf <- performance(preds, "tpr", "fpr")
plot(perf, colorize = TRUE)

#Считаем AUC
AUC.ROCR <- performance(preds, "auc")

# Добавляем прямую
abline(a = 0, b = 1)

# Добавляем значение AUC на график
text(x = .25, y = .65, paste("AUC = ", round(AUC.ROCR@y.values[[1]],5)))

# two.df = Subset(df_l, quality != 3)
# two.cv.qda <- qda(two.df[, 1:11], two.df[, 12], CV = T, prior=c(1,1)/2)
# predcv <- prediction(two.cv.qda$posterior[,2], two.df[, 12])
# perfcv <- performance(predcv, "tpr", "fpr")
# plot(perfcv, colorize = TRUE)
# AUCcv <- performance(predcv, "auc")
# abline(a = 0, b = 1)
# text(x = .25, y = .65, paste("AUC = ", round(AUCcv@y.values[[1]],5)))

Теперь возьмем класс \(2\) и \(3\).

# Подготовка данных
two.df = Subset(df_l, quality != 1)
two.qda <- qda(two.df[, 1:11], grouping = two.df[, 12], prior=c(1,1)/2)
two.qdap <- predict(two.qda, two.df[, 1:11])

# Переводит данные в стандартизированный формат. Первый параметр - предсказания, второй - настоящая классификация.
preds <- prediction(two.qdap$posterior[,2], two.df[, 12])

# Строим ROC-кривую и рисуем ее
perf <- performance(preds, "tpr", "fpr")
plot(perf, colorize = TRUE)

#Считаем AUC
AUC.ROCR <- performance(preds, "auc")

# Добавляем прямую
abline(a = 0, b = 1)

# Добавляем значение AUC на график
text(x = .25, y = .65, paste("AUC = ", round(AUC.ROCR@y.values[[1]],5)))

# two.df = Subset(df_l, quality != 1)
# two.cv.qda <- qda(two.df[, 1:11], two.df[, 12], CV = T, prior=c(1,1)/2)
# predcv <- prediction(two.cv.qda$posterior[,2], two.df[, 12])
# perfcv <- performance(predcv, "tpr", "fpr")
# plot(perfcv, colorize = TRUE)
# AUCcv <- performance(predcv, "auc")
# abline(a = 0, b = 1)
# text(x = .25, y = .65, paste("AUC = ", round(AUCcv@y.values[[1]],5)))

И последняя пара - классы \(1\) и \(3\).

# Подготовка данных
two.df = Subset(df_l, quality != 2)
two.qda <- qda(two.df[, 1:11], grouping = two.df[, 12], prior=c(1,1)/2)
two.qdap <- predict(two.qda, two.df[, 1:11])

# Переводит данные в стандартизированный формат. Первый параметр - предсказания, второй - настоящая классификация.
preds <- prediction(two.qdap$posterior[,2], two.df[, 12])

# Строим ROC-кривую и рисуем ее
perf <- performance(preds, "tpr", "fpr")
plot(perf, colorize = TRUE)

#Считаем AUC
AUC.ROCR <- performance(preds, "auc")

# Добавляем прямую
abline(a = 0, b = 1)

# Добавляем значение AUC на график
text(x = .25, y = .65, paste("AUC = ", round(AUC.ROCR@y.values[[1]],5)))

# two.df = Subset(df_l, quality != 2)
# two.cv.qda <- qda(two.df[, 1:11], two.df[, 12], CV = T, prior=c(1,1)/2)
# predcv <- prediction(two.cv.qda$posterior[,2], two.df[, 12])
# perfcv <- performance(predcv, "tpr", "fpr")
# plot(perfcv, colorize = TRUE)
# AUCcv <- performance(predcv, "auc")
# abline(a = 0, b = 1)
# text(x = .25, y = .65, paste("AUC = ", round(AUCcv@y.values[[1]],5)))